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Abstract. 

Debris disk images show clumps, rings, warps, and other structures, many of which have been 
interpreted as perturbations from hidden planets. But so far, no models of these structures have 
properly accounted for collisions between dust grains. We have developed new steady-state 3D 
models of debris disks that self-consistently incorporate grain-grain collisions. We summarize our 
algorithm and use it to illustrate how collisions interact with resonant trapping in the presence of a 
planet. 
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INTRODUCTION 

Debris disk images have emerged as important tools for finding extrasolar planets. A 
2005 image of Fomalhaut's debris disk [7] showed a ring offset from the central star 
pointing to a hidden planet; two years later, new HST Advanced Camera for Surveys 
(ACS) images revealed one of the first-ever directly imaged exoplanets orbiting in- 
side the disk |@]. The new generation of high contrast coronagraphs (NICI, HiCIAO, 
SPHERE, etc.) and mid-infrared/submillimeter imaging tools (JWST, SOFIA, ALMA) 
promises to return even more spectacular images of disks and planetary systems. But we 
cannot yet reliably decode the spectacular images of debris disks made by HST, JCMT, 
and other observatories, partly because no current models of structures in debris disks 
properly account for dust grain collisions. Many of the structures we see probably arise 
from the interaction between planetary perturbations, grain-grain, and grain-rock colli- 
sions. So even for signatures like the offset in the Fomalhaut ring that seem to be clear 
planet signposts, we have not yet been able to reliably derive planet orbital parameters. 

The key parameter that shows the relative importance of collisions for a population of 
grains in a debris disk is the ratio of the collision time, tcoii, to the Poynting-Robertson 
(PR) time, tpR, the time it takes the grains to spiral into the star under PR drag [2]. This 
ratio, rj = tpn/tcoiu represents roughly the number of collisions a particle undergoes 
during its lifetime in the absence of resonant trapping (which can increase the particle's 
lifetime by an order of magnitude.) The collision time, in turn, is approximately set by 
the cloud's face-on optical depth, T^; tcoii ~ torbu/^'^'^L where torbu is the orbital time 
[e.g. For typical imaged debrs disks (e.g. Vega, j8 Pictoris, Fomalhaut, etc.) r] is 
roughly 0.1-10 for a 10 jim grain. 

We describe here a new kind of debris disk model that self-consistently incorpo- 
rates grain-grain collisions. Our algorithm yields a dust distribution that simultaneously 
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FIGURE 1. Surface density as a function of rjo for a disk with an Earth-mass planet on a circular orbit 
at 1 AU, no fragmentation. The top row shows the entire disk, which extends to 4.25 AU. The bottom row 
shows zoomed-in views of the resonant ring. Collisions emphasize the birth ring and change the shape of 
the resonant ring. 



solves the governing dynamical equations and the mass flux equation in 3D. It uses the 
results of a coUisionless disk simulation as input and iteratively "grooms" the disk un- 
til it converges on the correct density distribution for a coUisional disk, allowing us to 
model disks with rj up to at least a few hundred. 



NUMERICAL TECHNIQUES 

To model the structures of dust clouds containing planets, we begin by constructing a 
"seed model". We numerically integrate the orbits of dust particles, a planet, and a star, 
as many authors have done 10, [121 191 M, H, [H, [13] • To the gravitational forces on the 
particle we add PR drag, radiation pressure, and solar wind drag. Radiation pressure 
effectively reduces the mass of the star as seen by the particle by the factor (1 — /3). 
In the solar system, /3 is roughly (0.285/im) /s for a spherical particle with radius s and 
density 2 g cm^-^. To model solar wind drag, we modify the PR drag coefficient, as many 
of the above authors have done. We release the particle in an orbit representing that of an 
asteroid or other parent body. The particle jumps to a new 2-body orbit appropriate for 
the radiation pressure it feels, conserving its velocity at release. We follow the particle 
until it moves outside our region of interest or collides with the Sun or a planet. At 
regular intervals, we record the location of the particle. We use the recorded locations 
to produce a 3D histogram modeling the steady-state distribution of particles released at 
regular intervals from a belt of small bodies. 



Our simulations use many more particles than the previous generation of models. 
This capability is important because the density distributions produced by the standard 
techniques suffer from three key sources of uncertainty that can be alleviated using 
larger numbers of particles: Poisson noise in each cell of the histogram, Poisson noise 
in the population of the mean motion resonances (MMRs) and the domination of the 
histograms by a few long-lived particles [e.g. 9]. Our models are the first to properly 
manage all three sources of noise. We handle Poisson noise by using at least 5,000 
particles per size bin, choosing an appropriate interval for recording particle positions 
and by devoting equal amounts of computer time to particles with large and small /3 
values. Our collision algorithms (discussed below) manage the third source of noise by 
destroying long-lived particles. 

To make the most of the available computers, we developed our own hybrid symplec- 
tic integrator adapted for the dust problem. This integrator combines a Burlisch-Stoer 
integrator and a symplectic integrator [13] in the manner of [3]; the Bulirsch-Stoer in- 
tegrator turns on during close encounters. We have been running our codes on NASA's 
Discover cluster (2,560 processors: 512 3.2 GHz plus 2,048 2.6 GHz). We maintain a 
free online database of our coUisionless debris disk simulations at 

http://asd.gsfc.nasa. gov/ Christopher . Stark /catalog. phpl 

Once the seed model is finished, we apply our new algorithm for incorporating grain- 
grain collisions. We use the following iterative approach: 

1. Integrate the orbits of the particles as described above (the seed model). 

2. Record both the positions and the velocities in a 6-dimensional histogram: a distri- 
bution function. 

3. Allow the particles to interact with the distribution function assembled in steps 1 
and 2; update the distribution function. 

4. Repeat step 3 until the process converges. 

Using the assumption that all collisions lead to complete, instantaneous grain removal 
we have produced the first 3D models of exozodiacal clouds that incorporate colli- 
sions in steady-state cloud self-consistently. 

In this picture, we view each of the "particles" in the simulation as particle streams 
representing billions of particles each. As they progress through the cloud, the streams 
become attenuated by collisions with other streams as e^^, where T is the optical 
depth of the cloud along the stream; T = / < no > ds. Between any two records in 
the histogram, we approximate the optical depth for collisions along the stream i as 
At ~ Y^jnj(yj\vi — Vj|A?, where rij, Oj and vj are the number density, particle cross 
section, and vector velocities of all the streams that flow through a given histogram 
cell, as measured in that cell, and At is the time interval between the two records. This 
approximation works as long At << tcoii, which is always true for our simulations. We 
use the iteration scheme described above to make sure that all the streams are attenuated 
consistently with one another. We examined five possible iterative schemes that achieve 
self-consistency. Two that met our requirements for accuracy; we chose the faster one. 
This scheme allows us to process a simulation containing 5,000 particle streams in 15 
minutes on a single processor. The algorithm typically converges in ~ 5 iterations to a 
point where the solution found in each subsequent iteration differs by no more than 1% 



from the previous iteration's solution in any cell. The time it takes to run a 5,000-particle 
simulation with our algorithm is dominated by the time it takes to integrate the orbits of 
the particles, not the collision algorithm. 

To test our simulations, we compared them to analytic solutions to the azimuthally- 
symmetric steady-state mass flux equation for dust grains under PR drag assuming that 
collisions serve only to completely remove particles [14]. We found that our simulations 
disagree substantially with the analytic solutions only for tjq > 235; for such short 
collision times, the histogram bins become too small to resolve the radial structure. 
To prevent this problem, we must simply be careful to choose the bin size, d, so that 
d<ro/{2rio). 

RESULTS AND DISCUSSION 

Figure 1 shows some examples of our first generation of coUisional disk models with no 
daughter particles. It shows the same disk/planet combination with three different values 
of rj. An Earth-mass planet in this simulation (semi-major axis 1 AU) traps some dust 
in MMRs, creating a resonant ring. An emerging common description of debris disks is 
that they have three radial regions ifllll 

I. A "birth ring" of large bodies, perhaps analogous to the asteroid belt or Kuiper 
belt. 

II. A region interior to the birth ring dominated by large grains spiraling in under 
PR drag and 

III. A skirt of small unbound or marginally -bound grains (/3-meteoroids) exterior 
to the birth ring. 

The simulation shown in Figure 1 models regions I and II. As rj increases, the birth 
ring becomes emphasized with respect to the resonant ring, which is in region II. The 
resonant ring also changes shape as rj increases; sharp structures in the inner part of the 
ring vanish, and the azimuthal structure shifts. Clearly, understanding these effects will 
be important for understanding images of debris disks and probably our own zodiacal 
cloud. Figure 2 shows the collision rate as a function of semi-major axis in one of our 
coUisional models (green). The collision rate is enhanced within the MMRs. The effect 
is even more pronounced in the disk images than this plot might suggest, because the 
particles spend most of their time in the MMRs. The first-order MMRs show double- 
peaked profiles, probably associated with their large libration amplitudes. An orange 
dashed curve shows the approximation, t^oii ~ torbit/^'^'^i used in previous work. 

So far, these results mostly illustrate the potential of this algorithm. We are working on 
new coUisional models that incorporate daughter particles, a range of particle sizes, and 
a range of parent body dynamics. We expect that these kinds of models will soon help 
answer some important questions in debris disk modeling: How robust are the structures 
predicted by 3-body simulations? How does the distribution of parent bodies impact the 
appearance of the dust cloud? They should be useful for recognizing planet migration 
in debris disks and providing better constraints on the orbital parameters of planets in 
eccentric rings like Fomalhaut. 
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FIGURE 2. Collision rate in one of our simulations vs. semimajor axis (green) compared to the naive 
estimate sometimes used in previous work (orange dashed). 



ACKNOWLEDGMENTS 

We thank NASA's Astrobiology program, NASA's Graduate Student Researchers Pro- 
gram and the NASA High-End Computing Program. 

REFERENCES 

1. Backman, D.E. & Paresce, F. 1993, In: Protostars and Planets III, p. 1253, Main-sequence stars with 
circumstellar solid material - The VEGA phenomenon 

2. Burns, J. A., Lamy, P. L., & Soter, S. 1979, Icarus, 40, 1, Radiation forces on small particles in the 
solar system 

3. Chambers, J.E. 1999, MNRAS, 304, 793, A Hybrid Symplectic Integrator That Permits Close Encoun- 
ters Between Massive Bodies 

4. Deller, A.T. & Maddison, S.T. 2005, ApJ, 625, 398, Numerical ModeHng of Dusty Debris Disks. 

5. Dermott, S.F. et al. 1994, Nature 369, 7 19, A circumsolar ring of asteroidal dust in resonant lock with 
the Earth. 

6. Kalas, P., et al. 2008, Science, 322, 1345, Optical Images of an Exosolar Planet 25 Light- Years from 
Earth 

7. Kalas, P., Graham, J. R., & Clampin, M. C, 2005, Nature, 435, 1067, A Planetary System as the Origin 
of Structure in FomalhautOs Dust Belt 

8. Moran, S.M., Kuchner, M.J, & Holman, M.J. 2004, ApJ, 612, 1163, The Dynamical Influence of a 
Planet at Semimajor Axis 3.4 AU on the Dust around ? Eridani 

9. Moro-Martin, A. & Malhotra, R. 2002, AJ, 124, 2305, A study of the dynamics of dust from the Kuiper 
belt: spatial distribution and spectral energy distribution. 

10. Stark, C. & Kuchner, M. J. 2008, ApJ, 686, 637, The DetectabiHty of Exo-Earths and Super-Earths 
Via Resonant Signatures in Exozodiacal Clouds 

11. Strubbe, L.E., & Chiang, E.I., 2006, ApJ, 648, 652, Dust Dynamics, Surface Brightness Profiles, and 
Thermal Spectra of Debris Disks: The Case of AU Microscopii 

12. Wilner, D.J., Holman, M.J, Kuchner, M.J, & Ho, PT.P 2002, ApJ, 569, LI 15-Ll 19, Structure in the 
Dusty Debris around Vega. 

13. Wisdom, J., & Holman, M. J. 1991, AJ, 102, 1528, Symplectic maps for the n-body problem 

14. Wyatt, M.C. 2005, A&A 433, 1007-1012, The insignificance of P-R drag in detectable extrasolar 
planetesimal belts. 

15. Wyatt, M.C. 2006, ApJ, 639, 1153, Dust in Resonant Extrasolar Kuiper Belts: Grain Size and 
Wavelength Dependence of Disk Structure 



